{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"./images/CC-BY.png\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli</td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "from __future__ import absolute_import\n",
    "\n",
    "%matplotlib inline\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Numerical Differentiation\n",
    "\n",
    "**GOAL:**  Given a set of $N+1$ points $(x_i, y_i)$ compute the derivative of a given order to a specified accuracy.\n",
    "\n",
    "**Approach:** Find the interpolating polynomial $P_N(x)$ and differentiate that."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Newton Polynomials\n",
    "\n",
    "While the interpolating polynomial $P_N(x)$ through $N+1$ distinct points is unique,  it can be expressed in terms of many bases.\n",
    "\n",
    "So far we have considered the monomial basis $< 1,x,x^2,\\ldots,x^n>$ and the Lagrange basis $\\ell_i(x)$\n",
    "\n",
    "However, for our purposes here it's somewhat more convenient to use the *Newton Polynomial* Basis\n",
    "\n",
    "\\begin{align}\n",
    "    n_0 &= 1\\\\\n",
    "    n_j(x) &= \\prod^{j-1}_{i=0} (x - x_i)\\quad\\mathrm{for}\\quad j>0 \\\\\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Newton Polynomials\n",
    "\n",
    "Unrolling them, we see\n",
    "\\begin{align}\n",
    "    n_0(x) &= 1\\\\\n",
    "    n_1(x) &= (x-x_0) \\\\\n",
    "    n_2(x) &= (x - x_0)(x-x_1)\\\\\n",
    "        &\\vdots \\\\\n",
    "    n_N(x) &= (x - x_0)(x - x_1)\\cdots(x-x_{N-1})\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Note:  The Newton polynomials have some features common to both the monomials and the Lagrange polynomials\n",
    "\n",
    "* Like the monomials they are monic polynomials of increasing degree in $x$, i.e. $n_2(x)$ is quadratic\n",
    "* Like the Lagrange polynomials, some Newton polynomials vanish identically at the interpolation points $x_i$.  In particular\n",
    "\n",
    "$$\n",
    "    n_j(x_i) = 0 \\quad \\mathrm{for}\\quad j> i\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Newton's Form\n",
    "\n",
    "This last feature makes it convenient to solve for the interpolating polynomial.\n",
    "\n",
    "$$P_N(x) = \\sum^N_{j=0} a_j n_j(x)$$\n",
    "\n",
    "as the $(N+1)\\times (N+1)$ linear system\n",
    "\n",
    "$$\n",
    "    P_N(x_i) = y_i\\quad\\mathrm{for}\\quad i=0,\\ldots,N\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Becomes the lower triangular system\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "1 & 0 & 0 & \\ldots & 0 \\\\\n",
    "1 & (x_1-x_0) & 0 & \\ldots & 0 \\\\\n",
    "1 & (x_2-x_0) & (x_2-x_0)(x_2-x_1) & \\ldots & 0 \\\\\n",
    "\\vdots & \\vdots & \\vdots & \\vdots & \\vdots  \\\\\n",
    "1 & n_1(x_n) & n_2(x_n) & \\ldots & n_n(x_n) \\\\\n",
    "\\end{bmatrix}\\begin{bmatrix}\n",
    "a_0 \\\\\n",
    "a_1 \\\\\n",
    "a_2 \\\\\n",
    "\\vdots \\\\\n",
    "a_n \\\\\n",
    "\\end{bmatrix}=\n",
    "\\begin{bmatrix}\n",
    "y_0 \\\\\n",
    "y_1 \\\\\n",
    "y_2 \\\\\n",
    "\\vdots \\\\\n",
    "y_n \\\\\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Which can be solved sequentially starting with\n",
    "\n",
    "\\begin{align}\n",
    "    a_0 & = y_0\\\\\n",
    "    a_1 & = \\frac{y_1 - y_0}{x_1 - x_0} \\\\\n",
    "    a_2 & = \\frac{y_2 - y_0 - a_1(x_2 - x_0)}{(x_2-x_0)(x_2-x_1)} \\\\\n",
    "    &\\vdots \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The divided differences\n",
    "\n",
    "The $a_j = [y_0, \\ldots, y_j]$ are also known as *the divided differences* which are defined recursively as\n",
    "\n",
    "$$[y_i] = y_i \\quad i \\in \\{0,\\ldots, N+1\\}$$\n",
    "\n",
    "and\n",
    "\n",
    "$$[y_i, \\ldots , y_{i+j}] = \\frac{[y_{i+1}, \\ldots , y_{i + j}] - [y_{i},\\ldots,y_{i+j-1}]}{x_{i+j} - x_{i}} \\quad i \\in \\{0,\\ldots,N+1 - j\\} \\quad j \\in \\{1,\\ldots, N+1\\}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "These formulas are recursively defined but not so helpful, here are a few examples to start out with:\n",
    "\n",
    "$$[y_0] = y_0$$\n",
    "\n",
    "$$[y_0, y_1] = \\frac{y_1 - y_0}{x_1 - x_0}$$\n",
    "\n",
    "$$[y_0, y_1, y_2] = \\frac{[y_1, y_2] - [y_0, y_1]}{x_{2} - x_{0}} = \\frac{\\frac{y_2 - y_1}{x_2 - x_1} - \\frac{y_1 - y_0}{x_1 - x_0}}{x_2 - x_0} = \\frac{y_2 - y_1}{(x_2 - x_1)(x_2 - x_0)} - \\frac{y_1 - y_0}{(x_1 - x_0)(x_2 - x_0)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The benefit of writing a polynomial like this is that it isolates the $x$ dependence (we can easily take derivatives of this form).\n",
    "\n",
    "In general then $P_N(x)$ can be written in Newton's form as\n",
    "\n",
    "$$P_N(x) = y_0 + (x-x_0)[y_0, y_1] + (x - x_0) (x - x_1) [y_0, y_1, y_2] + \\cdots + (x-x_0) (x-x_1) \\cdots (x-x_{N-1}) [y_0, y_1, \\ldots, y_{N}]$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def divided_difference(x, y, N=50):\n",
    "    \"\"\"Compute the Nth divided difference using *x* and *y*\"\"\"\n",
    "    if N == 0:\n",
    "        raise Exception(\"Reached recursion limit!\")\n",
    "    \n",
    "    # Reached the end of the recurssion\n",
    "    if y.shape[0] == 1:\n",
    "        return y[0]\n",
    "    elif y.shape[0] == 2:\n",
    "        return (y[1] - y[0]) / (x[1] - x[0])\n",
    "    else:\n",
    "        return (divided_difference(x[1:], y[1:], N=N-1) - divided_difference(x[:-1], y[:-1], N=N-1)) / (x[-1] - x[0])\n",
    "    \n",
    "def newton_basis(x, x0):\n",
    "    \"\"\" construct newton basis for the interpolating polynomial given interpolation points x0\"\"\"\n",
    "    basis = numpy.ones((len(x0), len(x)))\n",
    "    for j in range(len(x0)):\n",
    "        for i in range(j):\n",
    "            basis[j, :] *= (x - x0[i])\n",
    "    return basis\n",
    "\n",
    "def P_newton(x, data):\n",
    "    \"\"\" construct the interpolating polynomial of order N through the N+1 data points data using a Newton basis\"\"\"\n",
    "    P = numpy.zeros(x.shape)\n",
    "    basis = newton_basis(x, data[:,0])\n",
    "    for j in range(data.shape[0]):\n",
    "        P += divided_difference(data[:j + 1, 0], data[:j + 1, 1]) * basis[j, :]\n",
    "        \n",
    "    return P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Calculate a polynomial in Newton Form\n",
    "data = numpy.array([[-2.0, 1.0], [-1.5, -1.0], [-0.5, -3.0], [0.0, -2.0], [1.0, 3.0], [2.0, 1.0]])\n",
    "x = numpy.linspace(-2.0, 2.0, 100)\n",
    "basis = newton_basis(x, data[:,0])\n",
    "P = P_newton(x, data)\n",
    "N = data.shape[0] - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Plot basis and interpolant\n",
    "fig = plt.figure(figsize = (16,6))\n",
    "fig.subplots_adjust(hspace=.5)\n",
    "\n",
    "axes = [None, None]\n",
    "axes[0] = fig.add_subplot(1, 2, 1)\n",
    "axes[1] = fig.add_subplot(1, 2, 2)\n",
    "\n",
    "for j in range(N + 1):\n",
    "    axes[0].plot(x, basis[j, :])\n",
    "    axes[1].plot(data[j, 0], data[j, 1],'ko')\n",
    "axes[1].plot(x, P)\n",
    "\n",
    "axes[0].set_title(\"Newton Polynomial Basis\")\n",
    "axes[0].set_xlabel(\"x\")\n",
    "axes[0].set_ylabel(\"$n_j(x)$\")\n",
    "axes[0].grid()\n",
    "\n",
    "axes[1].set_title(\"Interpolant $P_%s(x)$\" % N)\n",
    "axes[1].set_xlabel(\"x\")\n",
    "axes[1].set_ylabel(\"$P_%s(x)$\" % N)\n",
    "axes[1].grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As another concrete example consider the general quadratic polynomial written in Newton's form\n",
    "\n",
    "$$P_2(x) = [y_0] + (x - x_0) [y_0, y_1] + (x - x_0)(x - x_1) [y_0, y_1, y_2] \\\\= y_0 + (x - x_0) \\frac{y_1 - y_0}{x_1 - x_0} + (x - x_0)(x - x_1) \\left ( \\frac{y_2 - y_1}{(x_2 - x_1)(x_2 - x_0)} - \\frac{y_1 - y_0}{(x_1 - x_0)(x_2 - x_0)} \\right )$$\n",
    "\n",
    "Recall that the interpolating polynomial of degree $N$ through these points is unique!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Error Analysis\n",
    "\n",
    "Given $N + 1$ points we can form an interpolant $P_N(x)$ of degree $N$ where\n",
    "\n",
    "$$f(x) = P_N(x) + R_N(x)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "And we know from Lagrange's Theorem that the remainder term looks like\n",
    "\n",
    "$$R_N(x) = (x - x_0)(x - x_1)\\cdots (x - x_{N})(x - x_{N+1}) \\frac{f^{(N+1)}(c)}{(N+1)!}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "noting that we need to require that $f(x) \\in C^{N+1}$ on the interval of interest.  Taking the derivative of the interpolant $P_N(x)$ then leads to \n",
    "\n",
    "$$P_N'(x) = [y_0, y_1] + ((x - x_1) + (x - x_0)) [y_0, y_1, y_2] + \\cdots + \\left(\\sum^{N-1}_{i=0}\\left( \\prod^{N-1}_{j=0,~j\\neq i} (x - x_j) \\right )\\right ) [y_0, y_1, \\ldots, y_N]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Similarly we can find the derivative of the remainder term $R_N(x)$ as\n",
    "\n",
    "$$R_N'(x) = \\left(\\sum^{N}_{i=0} \\left( \\prod^{N}_{j=0,~j\\neq i} (x - x_j) \\right )\\right ) \\frac{f^{(N+1)}(c)}{(N+1)!}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now if we consider the approximation of the derivative evaluated at one of our data points $(x_k, y_k)$ these expressions simplify such that\n",
    "\n",
    "$$f'(x_k) = P_N'(x_k) + R_N'(x_k)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "If we let $\\Delta x = \\max_i |x_k - x_i|$ we then know that the remainder term will be $\\mathcal{O}(\\Delta x^N)$ as $\\Delta x \\rightarrow 0$ thus showing that this approach converges and we can find arbitrarily high order approximations (ignoring floating point error)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### <font color='red'>Caution</font>\n",
    "\n",
    "High order does not necessarily imply high-accuracy!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def PN_prime(x, data):\n",
    "    \"\"\" calculate derivative of interpolating polynomail at point x\"\"\"\n",
    "    # General form of derivative of P_N'(x)\n",
    "    P_prime = numpy.zeros(x.shape)\n",
    "    newton_basis_prime = numpy.empty(x.shape)\n",
    "    product = numpy.empty(x.shape)\n",
    "    N = data.shape[0] - 1\n",
    "    for n in range(N + 1):\n",
    "        newton_basis_prime = 0.0\n",
    "        for i in range(n):\n",
    "            product = 1.0\n",
    "            for j in range(n):\n",
    "                if j != i:\n",
    "                    product *= (x - data[j, 0])\n",
    "            newton_basis_prime += product\n",
    "        P_prime += divided_difference(data[:n+1, 0], data[:n+1, 1]) * newton_basis_prime\n",
    "        \n",
    "    return P_prime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "num_points = 4\n",
    "xmin = -2.\n",
    "xmax = 2.\n",
    "data = numpy.empty((num_points, 2))\n",
    "data[:, 0] = numpy.linspace(xmin, xmax, num_points)\n",
    "data[:, 1] = numpy.sin(data[:, 0])\n",
    "x = numpy.linspace(xmin, xmax, 100)\n",
    "P_prime = PN_prime(x, data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, numpy.cos(x), 'k', label=\"$f'(x)$\")\n",
    "axes.plot(x, P_prime, 'r--', label=\"$Pn'(x)$\")\n",
    "axes.plot(data[:,0], numpy.cos(data[:,0]), 'ro')\n",
    "axes.legend(loc='best')\n",
    "axes.set_title(\"$f'(x)$\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$f'(x)$ and $\\hat{f}'(x)$\")\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Examples\n",
    "\n",
    "Often in practice we only use a small number of data points to derive a differentiation formula.  In the context of differential equations we also often have $f(x)$ so that $f(x_k) = y_k$ and we can approximate the derivative of a known function $f(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example 1:  1st order Forward and Backward Differences\n",
    "\n",
    "Using 2 points we can get an approximation that is $\\mathcal{O}(\\Delta x)$:\n",
    "\n",
    "$$f'(x) \\approx P_1'(x) = [y_0, y_1] = \\frac{y_1 - y_0}{x_1 - x_0} = \\frac{y_1 - y_0}{\\Delta x} = \\frac{f(x_1) - f(x_0)}{\\Delta x}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can also calculate the error as\n",
    "\n",
    "$$R_1'(x) = -\\Delta x \\frac{f''(c)}{2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can also derive the \"forward\" and \"backward\" formulas by considering the question slightly differently.  Say we want $f'(x_n)$, then the \"forward\" finite-difference can be written as\n",
    "\n",
    "$$f'(x_n) \\approx D_1^+ = \\frac{f(x_{n+1}) - f(x_n)}{\\Delta x}$$\n",
    "\n",
    "and the \"backward\" finite-difference as\n",
    "\n",
    "$$f'(x_n) \\approx D_1^- = \\frac{f(x_n) - f(x_{n-1})}{\\Delta x}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Note these approximations should be familiar to use as the limit as $\\Delta x \\rightarrow 0$ these are no longer approximations but equivalent definitions of the derivative at $x_n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: numpy.sin(x)\n",
    "f_prime = lambda x: numpy.cos(x)\n",
    "\n",
    "# Use uniform discretization\n",
    "x = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, 1000)\n",
    "N = 20\n",
    "x_hat = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, N)\n",
    "delta_x = x_hat[1] - x_hat[0]\n",
    "\n",
    "# Compute forward difference using a loop\n",
    "f_prime_hat = numpy.empty(x_hat.shape)\n",
    "for i in range(N - 1):\n",
    "    f_prime_hat[i] = (f(x_hat[i+1]) - f(x_hat[i])) / delta_x\n",
    "f_prime_hat[-1] = (f(x_hat[-1]) - f(x_hat[-2])) / delta_x\n",
    "\n",
    "# Vector based calculation\n",
    "# f_prime_hat[:-1] = (f(x_hat[1:]) - f(x_hat[:-1])) / (delta_x)\n",
    "\n",
    "# Use first-order differences for points at edge of domain\n",
    "f_prime_hat[-1] = (f(x_hat[-1]) - f(x_hat[-2])) / delta_x  # Backward Difference at x_N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(x, f_prime(x), 'k')\n",
    "axes.plot(x_hat + 0.5 * delta_x, f_prime_hat, 'ro')\n",
    "axes.set_xlim((x[0], x[-1]))\n",
    "axes.set_ylim((-1.1, 1.1))\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Aside:  Computing Order of Convergence\n",
    "\n",
    "Say we had the error $E(\\Delta x)$ and we wanted to make a statement about the rate of convergence (note we can replace $E$ here with the $R$ from above).  Then we can do the following:\n",
    "$$\\begin{aligned}\n",
    "    E(\\Delta x) &= C \\Delta x^n \\\\\n",
    "    \\log E(\\Delta x) &= \\log C + n \\log \\Delta x\n",
    "\\end{aligned}$$\n",
    "\n",
    "The slope of the line is $n$ when modeling the error like this!  We can also match the first point by solving for $C$:\n",
    "\n",
    "$$\n",
    "    C = e^{\\log E(\\Delta x) - n \\log \\Delta x}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: numpy.sin(x)\n",
    "f_prime = lambda x: numpy.cos(x)\n",
    "\n",
    "# Compute the error as a function of delta_x\n",
    "N_range = numpy.logspace(1, 4, 10, dtype=int)\n",
    "delta_x = numpy.empty(N_range.shape)\n",
    "error = numpy.empty((N_range.shape[0], 4))\n",
    "for (i, N) in enumerate(N_range):\n",
    "    x_hat = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, N)\n",
    "    delta_x[i] = x_hat[1] - x_hat[0]\n",
    "\n",
    "    # Compute forward difference\n",
    "    f_prime_hat = numpy.empty(x_hat.shape)\n",
    "    f_prime_hat[:-1] = (f(x_hat[1:]) - f(x_hat[:-1])) / (delta_x[i])\n",
    "\n",
    "    # Use first-order differences for points at edge of domain\n",
    "    f_prime_hat[-1] = (f(x_hat[-1]) - f(x_hat[-2])) / delta_x[i]  # Backward Difference at x_N\n",
    "    \n",
    "    # The differences in error computations is interesting here.  Note that the L_\\infty norm returns a single\n",
    "    # point-wise like error where as the L_2 and L_1 are adding error up so they exhbit less convergence\n",
    "    error[i, 0] = numpy.linalg.norm(numpy.abs(f_prime(x_hat) - f_prime_hat), ord=numpy.infty)\n",
    "    error[i, 1] = numpy.linalg.norm(numpy.abs(f_prime(x_hat + 0.5 * delta_x[i]) - f_prime_hat), ord=numpy.infty)\n",
    "    error[i, 2] = numpy.linalg.norm(numpy.abs(f_prime(x_hat) - f_prime_hat), ord=1)\n",
    "    error[i, 3] = numpy.linalg.norm(numpy.abs(f_prime(x_hat) - f_prime_hat), ord=2)\n",
    "    \n",
    "error = numpy.array(error)\n",
    "delta_x = numpy.array(delta_x)\n",
    "    \n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "fig.set_figheight(fig.get_figheight() * 2)\n",
    "# plt.rc('legend', fontsize=6)\n",
    "# plt.rc('font', size=6)\n",
    "\n",
    "error_type = ['$L^\\infty$', 'offset $L^\\infty$', '$L^1$', '$L^2$']\n",
    "for i in range(error.shape[1]):\n",
    "    axes = fig.add_subplot(2, 2, i + 1)\n",
    "    axes.loglog(delta_x, error[:, i], 'ko', label=\"Error\")\n",
    "\n",
    "    order_C = lambda delta_x, error, order: numpy.exp(numpy.log(error) - order * numpy.log(delta_x))\n",
    "    axes.loglog(delta_x, order_C(delta_x[0], error[0, i], 1.0) * delta_x**1.0, 'r--', label=\"1st Order\")\n",
    "    axes.loglog(delta_x, order_C(delta_x[0], error[0, i], 2.0) * delta_x**2.0, 'b--', label=\"2nd Order\")\n",
    "    axes.legend(loc=4)\n",
    "    axes.set_title(\"Convergence of 1st Order Differences - %s\" % error_type[i])\n",
    "    axes.set_xlabel(\"$\\Delta x$\")\n",
    "    axes.set_ylabel(\"$|f'(x) - \\hat{f}'(x)|$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example 2: 2nd Order Centered Difference\n",
    "\n",
    "Now lets use 3 points to calculate the 2nd order accurate finite-difference.  Consider the points $(x_{n}, y_{n})$, $(x_{n-1}, y_{n-1})$, and $(x_{n+1}, y_{n+1})$, from before we have\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    P_2(x) &= [f(x_0)] + (x - x_0) [f(x_0), f(x_1)] + (x - x_0)(x - x_1) [f(x_0), f(x_1), f(x_2)] \\\\\n",
    "    &= f(x_0) + (x - x_0) \\frac{f(x_1) - f(x_0)}{x_1 - x_0} + (x - x_0)(x - x_1) \\left ( \\frac{f(x_2) - f(x_1)}{(x_2 - x_1)(x_2 - x_0)} - \\frac{f(x_1) - f(x_0)}{(x_1 - x_0)(x_2 - x_0)} \\right )\n",
    "\\end{aligned}$$\n",
    "Compute the formula for the derivative.  Assume that the distance between the $x_i$ are equal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\\begin{aligned}\n",
    "    P_2'(x) &= [f(x_n), f(x_{n+1})] + ((x - x_n) + (x - x_{n+1})) [f(x_n), f(x_{n+1}), f(x_{n-1})] \\\\\n",
    "    &= \\frac{f(x_{n+1}) - f(x_n)}{x_{n+1} - x_n}  + ((x - x_n) + (x - x_{n+1}))\\\\ &* \\left ( \\frac{f(x_{n-1}) - f(x_{n+1})}{(x_{n-1} - x_{n+1})(x_{n-1} - x_n)} - \\frac{f(x_{n+1}) - f(x_n)}{(x_{n+1} - x_n)(x_{n-1} - x_n)} \\right )\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Evaluating at $x_n$ and assuming the points $x_{n-1}, x_n, x_{n+1}$ are evenly spaced leads to\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    P_2'(x_n) &= \\frac{f(x_{n+1}) - f(x_n)}{\\Delta x} - \\Delta x \\left ( \\frac{f(x_{n-1}) - f(x_{n+1})}{2\\Delta x^2} + \\frac{f(x_{n+1}) - f(x_n)}{\\Delta x^2} \\right ) \\\\\n",
    "    &=\\frac{f(x_{n+1}) - f(x_n)}{\\Delta x} - \\left ( \\frac{f(x_{n+1}) - 2f(x_n) + f(x_{n-1})}{2\\Delta x}\\right ) \\\\\n",
    "    &=\\frac{2f(x_{n+1}) - 2f(x_n) - f(x_{n+1}) + 2f(x_n) - f(x_{n-1})}{2 \\Delta x} \\\\\n",
    "    &=\\frac{f(x_{n+1}) - f(x_{n-1})}{2 \\Delta x}\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This finite-difference is second order accurate and is centered about the point it is meant to approximate ($x_n$).  We can show that it is second order by again considering the remainder term's derivative\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    R_2'(x) &= \\left(\\sum^{2}_{i=0} \\left( \\prod^{2}_{j=0,~j\\neq i} (x - x_j) \\right )\\right ) \\frac{f'''(c)}{3!} \\\\\n",
    "    &= \\left ( (x - x_{n+1}) (x - x_{n-1}) + (x-x_n) (x-x_{n-1}) + (x-x_n)(x-x_{n+1}) \\right ) \\frac{f'''(c)}{3!}\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Again evaluating this expression at $x = x_n$ and assuming evenly space points we have\n",
    "\n",
    "$$R_2'(x_n) = -\\Delta x^2 \\frac{f'''(c)}{3!}$$\n",
    "\n",
    "showing that our error is $\\mathcal{O}(\\Delta x^2)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: numpy.sin(x)\n",
    "f_prime = lambda x: numpy.cos(x)\n",
    "\n",
    "# Use uniform discretization\n",
    "x = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, 1000)\n",
    "N = 21\n",
    "x_hat = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, N)\n",
    "delta_x = x_hat[1] - x_hat[0]\n",
    "\n",
    "# Compute derivative\n",
    "f_prime_hat = numpy.empty(x_hat.shape)\n",
    "f_prime_hat[1:-1] = (f(x_hat[2:]) - f(x_hat[:-2])) / (2. * delta_x)\n",
    "\n",
    "# Use first-order differences for points at edge of domain\n",
    "f_prime_hat[0] = (f(x_hat[1]) - f(x_hat[0])) / delta_x     # Forward Difference at x_0\n",
    "f_prime_hat[-1] = (f(x_hat[-1]) - f(x_hat[-2])) / delta_x  # Backward Difference at x_N\n",
    "#f_prime_hat[0] = (-3.0 * f(x_hat[0]) + 4.0 * f(x_hat[1]) - f(x_hat[2])) / (2.0 * delta_x)\n",
    "#f_prime_hat[-1] = (3.0 * f(x_hat[-1]) - 4.0 * f(x_hat[-2]) + f(x_hat[-3])) / (2.0 * delta_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "plt.rcdefaults()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(x, f_prime(x), 'k',label=\"$f'(x)$\")\n",
    "axes.plot(x_hat, f_prime_hat, 'ro')\n",
    "axes.plot(x, f(x),'g--',label=\"$f(x)$\")\n",
    "axes.plot(x_hat, f(x_hat), 'go')\n",
    "axes.set_xlim((x[0], x[-1]))\n",
    "# axes.set_ylim((-1.1, 1.1))\n",
    "axes.grid()\n",
    "axes.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Compute the error as a function of delta_x\n",
    "delta_x = []\n",
    "error = []\n",
    "# for N in range(2, 101):\n",
    "for N in range(50, 1000, 50):\n",
    "    x_hat = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, N + 1)\n",
    "    delta_x.append(x_hat[1] - x_hat[0])\n",
    "\n",
    "    # Compute derivative\n",
    "    f_prime_hat = numpy.empty(x_hat.shape)\n",
    "    f_prime_hat[1:-1] = (f(x_hat[2:]) - f(x_hat[:-2])) / (2 * delta_x[-1])\n",
    "\n",
    "    # Use first-order differences for points at edge of domain\n",
    "    f_prime_hat[0] = (f(x_hat[1]) - f(x_hat[0])) / delta_x[-1]  \n",
    "    f_prime_hat[-1] = (f(x_hat[-1]) - f(x_hat[-2])) / delta_x[-1]\n",
    "    # Use second-order differences for points at edge of domain\n",
    "    f_prime_hat[0] =  (-3.0 * f(x_hat[0])  +  4.0 * f(x_hat[1])  + - f(x_hat[2]))  / (2.0 * delta_x[-1])\n",
    "    f_prime_hat[-1] = ( 3.0 * f(x_hat[-1]) + -4.0 * f(x_hat[-2]) +   f(x_hat[-3])) / (2.0 * delta_x[-1])\n",
    "    \n",
    "    # Compute Error\n",
    "    error.append(numpy.linalg.norm(numpy.abs(f_prime(x_hat) - f_prime_hat), ord=numpy.infty))\n",
    "#     error.append(numpy.linalg.norm(numpy.abs(f_prime(x_hat) - f_prime_hat), ord=2))\n",
    "    \n",
    "error = numpy.array(error)\n",
    "delta_x = numpy.array(delta_x)\n",
    "    \n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.loglog(delta_x, error, \"ro\", label=\"Approx. Derivative\")\n",
    "\n",
    "order_C = lambda delta_x, error, order: numpy.exp(numpy.log(error) - order * numpy.log(delta_x))\n",
    "axes.loglog(delta_x, order_C(delta_x[0], error[0], 1.0) * delta_x**1.0, 'b--', label=\"1st Order\")\n",
    "axes.loglog(delta_x, order_C(delta_x[0], error[0], 2.0) * delta_x**2.0, 'r--', label=\"2nd Order\")\n",
    "axes.legend(loc=4)\n",
    "axes.set_title(\"Convergence of 2nd Order Differences\")\n",
    "axes.set_xlabel(\"$\\Delta x$\")\n",
    "axes.set_ylabel(\"$|f'(x) - \\hat{f}'(x)|$\")\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example 3: Alternative Derivations\n",
    "\n",
    "An alternative method for finding finite-difference formulas is by using Taylor series expansions about the point we want to approximate.  The Taylor series about $x_n$ is\n",
    "\n",
    "$$f(x) = f(x_n) + (x - x_n) f'(x_n) + \\frac{(x - x_n)^2}{2!} f''(x_n) + \\frac{(x - x_n)^3}{3!} f'''(x_n) + \\mathcal{O}((x - x_n)^4)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Say we want to derive the second order accurate, first derivative approximation that we just did, this requires the values $(x_{n+1}, f(x_{n+1})$ and $(x_{n-1}, f(x_{n-1})$.  We can express these values via our Taylor series approximation above as\n",
    "\n",
    "\\begin{aligned}\n",
    "    f(x_{n+1}) &= f(x_n) + (x_{n+1} - x_n) f'(x_n) + \\frac{(x_{n+1} - x_n)^2}{2!} f''(x_n) + \\frac{(x_{n+1} - x_n)^3}{3!} f'''(x_n) + \\mathcal{O}((x_{n+1} - x_n)^4) \\\\\n",
    "\\end{aligned}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or\n",
    "\\begin{aligned}\n",
    "&= f(x_n) + \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) + \\frac{\\Delta x^3}{3!} f'''(x_n) + \\mathcal{O}(\\Delta x^4)\n",
    "\\end{aligned}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "and\n",
    "\n",
    "\\begin{align}\n",
    "f(x_{n-1}) &= f(x_n) + (x_{n-1} - x_n) f'(x_n) + \\frac{(x_{n-1} - x_n)^2}{2!} f''(x_n) + \\frac{(x_{n-1} - x_n)^3}{3!} f'''(x_n) + \\mathcal{O}((x_{n-1} - x_n)^4) \n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\\begin{align} \n",
    "&= f(x_n) - \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) - \\frac{\\Delta x^3}{3!} f'''(x_n) + \\mathcal{O}(\\Delta x^4)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Or all together (for regularly spaced points),\n",
    "\\begin{align} \n",
    "f(x_{n+1}) &= f(x_n) + \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) + \\frac{\\Delta x^3}{3!} f'''(x_n) + \\mathcal{O}(\\Delta x^4)\\\\\n",
    "f(x_{n-1})&= f(x_n) - \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) - \\frac{\\Delta x^3}{3!} f'''(x_n) + \\mathcal{O}(\\Delta x^4)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now to find out how to combine these into an expression for the derivative we assume our approximation looks like\n",
    "\n",
    "$$\n",
    "    f'(x_n) + R(x_n) = A f(x_{n+1}) + B f(x_n) + C f(x_{n-1})\n",
    "$$\n",
    "\n",
    "where $R(x_n)$ is our error.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Plugging in the Taylor series approximations we find\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    f'(x_n) + R(x_n) &= A \\left ( f(x_n) + \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) + \\frac{\\Delta x^3}{3!} f'''(x_n) + \\mathcal{O}(\\Delta x^4)\\right ) \\\\\n",
    "    & + B f(x_n) \\\\ \n",
    "    & + C \\left ( f(x_n) - \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) - \\frac{\\Delta x^3}{3!} f'''(x_n) + \\mathcal{O}(\\Delta x^4) \\right )\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Since we want $R(x_n) = \\mathcal{O}(\\Delta x^2)$ we want all terms lower than this to cancel except for those multiplying $f'(x_n)$ as those should sum to 1 to give us our approximation.  Collecting the terms with common evaluations of the derivatives on $f(x_n)$ we get a series of expressions for the coefficients $A$, $B$, and $C$ based on the fact we want an approximation to $f'(x_n)$.  The $n=0$ terms collected are $A + B + C$ and are set to 0 as we want the $f(x_n)$ term to also cancel.\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    f(x_n):&  &A + B + C &= 0 \\\\\n",
    "    f'(x_n): & &A \\Delta x - C \\Delta x &= 1 \\\\\n",
    "    f''(x_n): & &A \\frac{\\Delta x^2}{2} + C \\frac{\\Delta x^2}{2} &= 0\n",
    "\\end{aligned} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This last equation $\\Rightarrow A = -C$, using this in the second equation gives $A = \\frac{1}{2 \\Delta x}$ and $C = -\\frac{1}{2 \\Delta x}$.  The first equation then leads to $B = 0$.  Putting this altogether then gives us our previous expression including an estimate for the error:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    f'(x_n) + R(x_n) &= \\quad \\frac{1}{2 \\Delta x} \\left ( f(x_n) + \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) + \\frac{\\Delta x^3}{3!} f'''(x_n) + \\mathcal{O}(\\Delta x^4)\\right ) \\\\\n",
    "    & \\quad + 0 \\cdot f(x_n) \\\\ \n",
    "    & \\quad - \\frac{1}{2 \\Delta x} \\left ( f(x_n) - \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) - \\frac{\\Delta x^3}{3!} f'''(x_n) + \\mathcal{O}(\\Delta x^4) \\right ) \\\\\n",
    "    &=  f'(x_n) + \\frac{1}{2 \\Delta x} \\left ( \\frac{2 \\Delta x^3}{3!} f'''(x_n) + \\mathcal{O}(\\Delta x^4)\\right )\n",
    "\\end{aligned}$$\n",
    "so that we find\n",
    "$$\n",
    "    R(x_n) = \\frac{\\Delta x^2}{3!} f'''(x_n) + \\mathcal{O}(\\Delta x^3) = \\mathcal{O}(\\Delta x^2)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "#### Another way...\n",
    "\n",
    "There is one more way to derive the second order accurate, first order finite-difference formula.  Consider the two first order forward and backward finite-differences averaged together:\n",
    "\n",
    "$$\\frac{D_1^+(f(x_n)) + D_1^-(f(x_n))}{2} = \\frac{f(x_{n+1}) - f(x_n) + f(x_n) - f(x_{n-1})}{2 \\Delta x} = \\frac{f(x_{n+1}) - f(x_{n-1})}{2 \\Delta x}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example 4: Higher Order Derivatives\n",
    "\n",
    "Using our Taylor series approach lets derive the second order accurate second derivative formula.  Again we will use the same points and the Taylor series centered at $x = x_n$ so we end up with the same expression as before:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    f''(x_n) + R(x_n) &= \\quad A \\left ( f(x_n) + \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) + \\frac{\\Delta x^3}{3!} f'''(x_n) + \\frac{\\Delta x^4}{4!} f^{(4)}(x_n) + \\mathcal{O}(\\Delta x^5)\\right ) \\\\\n",
    "    &+ \\quad B \\cdot f(x_n) \\\\\n",
    "    &+ \\quad C \\left ( f(x_n) - \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) - \\frac{\\Delta x^3}{3!} f'''(x_n) + \\frac{\\Delta x^4}{4!} f^{(4)}(x_n) + \\mathcal{O}(\\Delta x^5) \\right )\n",
    "\\end{aligned}$$\n",
    "\n",
    "except this time we want to leave $f''(x_n)$ on the right hand side.  \n",
    "\n",
    "Try out the same trick as before and see if you can setup the equations that need to be solved."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Doing the same trick as before we have the following expressions:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    f(x_n): & & A + B + C &= 0\\\\\n",
    "    f'(x_n): & & A \\Delta x - C \\Delta x &= 0\\\\\n",
    "    f''(x_n): & & A \\frac{\\Delta x^2}{2} + C \\frac{\\Delta x^2}{2} &= 1\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The second equation implies $A = C$ which combined with the third implies\n",
    "\n",
    "$$A = C = \\frac{1}{\\Delta x^2}$$\n",
    "\n",
    "Finally the first equation gives\n",
    "\n",
    "$$B = -\\frac{2}{\\Delta x^2}$$\n",
    "\n",
    "leading to the final expression\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    f''(x_n) + R(x_n) &= \\quad \\frac{1}{\\Delta x^2} \\left ( f(x_n) + \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) + \\frac{\\Delta x^3}{3!} f'''(x_n) + \\frac{\\Delta x^4}{4!} f^{(4)}(x_n) + \\mathcal{O}(\\Delta x^5)\\right ) \\\\\n",
    "    &+ \\quad -\\frac{2}{\\Delta x^2} \\cdot f(x_n) \\\\\n",
    "    &+ \\quad \\frac{1}{\\Delta x^2} \\left ( f(x_n) - \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) - \\frac{\\Delta x^3}{3!} f'''(x_n) + \\frac{\\Delta x^4}{4!} f^{(4)}(x_n) + \\mathcal{O}(\\Delta x^5) \\right ) \\\\\n",
    "    &= f''(x_n) + \\frac{1}{\\Delta x^2} \\left(\\frac{2 \\Delta x^4}{4!} f^{(4)}(x_n) + \\mathcal{O}(\\Delta x^5) \\right )\n",
    "\\end{aligned}\n",
    "$$\n",
    "so that\n",
    "\n",
    "$$\n",
    "    R(x_n) = \\frac{\\Delta x^2}{12} f^{(4)}(x_n) + \\mathcal{O}(\\Delta x^3)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: numpy.sin(x)\n",
    "f_dubl_prime = lambda x: -numpy.sin(x)\n",
    "\n",
    "# Use uniform discretization\n",
    "x = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, 1000)\n",
    "N = 40\n",
    "x_hat = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, N)\n",
    "delta_x = x_hat[1] - x_hat[0]\n",
    "\n",
    "\n",
    "\n",
    "# Compute derivative\n",
    "f_dubl_prime_hat = numpy.empty(x_hat.shape)\n",
    "f_dubl_prime_hat[1:-1] = (f(x_hat[2:]) -2.0 * f(x_hat[1:-1]) + f(x_hat[:-2])) / (delta_x**2)\n",
    "\n",
    "# Use first-order differences for points at edge of domain\n",
    "f_dubl_prime_hat[0] = (2.0 * f(x_hat[0]) - 5.0 * f(x_hat[1]) + 4.0 * f(x_hat[2]) - f(x_hat[3])) / delta_x**2\n",
    "f_dubl_prime_hat[-1] = (2.0 * f(x_hat[-1]) - 5.0 * f(x_hat[-2]) + 4.0 * f(x_hat[-3]) - f(x_hat[-4])) / delta_x**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(x, f_dubl_prime(x), 'k')\n",
    "axes.plot(x_hat, f_dubl_prime_hat, 'ro')\n",
    "axes.set_xlim((x[0], x[-1]))\n",
    "axes.set_ylim((-1.1, 1.1))\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: numpy.sin(x)\n",
    "f_dubl_prime = lambda x: -numpy.sin(x)\n",
    "\n",
    "# Compute the error as a function of delta_x\n",
    "delta_x = []\n",
    "error = []\n",
    "# for N in xrange(2, 101):\n",
    "for N in range(50, 1000, 50):\n",
    "    x_hat = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, N)\n",
    "    delta_x.append(x_hat[1] - x_hat[0])\n",
    "\n",
    "    # Compute derivative\n",
    "    f_dubl_prime_hat = numpy.empty(x_hat.shape)\n",
    "    f_dubl_prime_hat[1:-1] = (f(x_hat[2:]) -2.0 * f(x_hat[1:-1]) + f(x_hat[:-2])) / (delta_x[-1]**2)\n",
    "\n",
    "    # Use second-order differences for points at edge of domain\n",
    "    f_dubl_prime_hat[0] = (2.0 * f(x_hat[0]) - 5.0 * f(x_hat[1]) + 4.0 * f(x_hat[2]) - f(x_hat[3])) / delta_x[-1]**2\n",
    "    f_dubl_prime_hat[-1] = (2.0 * f(x_hat[-1]) - 5.0 * f(x_hat[-2]) + 4.0 * f(x_hat[-3]) - f(x_hat[-4])) / delta_x[-1]**2\n",
    "    \n",
    "    error.append(numpy.linalg.norm(numpy.abs(f_dubl_prime(x_hat) - f_dubl_prime_hat), ord=numpy.infty))\n",
    "    \n",
    "error = numpy.array(error)\n",
    "delta_x = numpy.array(delta_x)\n",
    "    \n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "# axes.plot(delta_x, error)\n",
    "axes.loglog(delta_x, error, \"ko\", label=\"Approx. Derivative\")\n",
    "order_C = lambda delta_x, error, order: numpy.exp(numpy.log(error) - order * numpy.log(delta_x))\n",
    "axes.loglog(delta_x, order_C(delta_x[2], error[2], 1.0) * delta_x**1.0, 'b--', label=\"1st Order\")\n",
    "axes.loglog(delta_x, order_C(delta_x[2], error[2], 2.0) * delta_x**2.0, 'r--', label=\"2nd Order\")\n",
    "axes.legend(loc=4)\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
